<!DOCTYPE html PUBLIC "-//W3C//DTD HTML 4.01 Transitional//EN" "http://www.w3.org/TR/1999/REC-html401-19991224/loose.dtd">
<html lang="en">
<head>
	<meta http-equiv="Content-Type" content="text/html; charset=utf-8">
	<title>Attenuation Compensation Using Time Reversal Example (k-Wave)</title>
	<link rel="stylesheet" href="kwavehelpstyle.css" type="text/css">
	<meta name="description" content="Attenuation Compensation Using Time Reversal Example.">
</head>

<body><div class="content">

<h1>Attenuation Compensation Using Time Reversal Example</h1>

<p>This example demonstrates how the acoustic attenuation present in the photoacoustic forward problem can be compensated for using time reversal image reconstruction. It builds on the <a href="example_pr_2D_tr_circular_sensor.html">2D Time Reversal Reconstruction For A Circular Sensor Example</a>.</p>

<p>
    <ul>
        <li><a href="matlab:edit([getkWavePath('examples') 'example_pr_2D_TR_absorption_compensation.m']);" target="_top">Open the file in the MATLAB Editor</a></li>
        <li><a href="matlab:run([getkWavePath('examples') 'example_pr_2D_TR_absorption_compensation']);" target="_top">Run the file in MATLAB</a></li>
    </ul>
</p>

<p>For a more detailed discussion of this example and the underlying techniques, see B. E. Treeby, E. Z. Zhang, and B. T. Cox, "Photoacoustic tomography in absorbing acoustic media using time reversal," <em>Inverse Problems</em>, vol. 26, no. 11, p. 115003, 2010.</p>

<h2>Contents</h2>
<div>
	<ul>
		<li><a href="#heading2">Running the forward simulation</a></li>
		<li><a href="#heading3">Selecting the regularisation parameters</a></li>
		<li><a href="#heading4">Time reversal image reconstruction</a></li>
	</ul>
</div>

<a name="heading2"></a>
<h2>Running the forward simulation</h2>

<p>The sensor data is simulated using <code><a href="kspaceFirstOrder2D.html">kspaceFirstOrder2D</a></code> in a similar manner to previous examples. The initial pressure is set to the Shepp Logan phantom and a circular Cartesian sensor mask with 200 sensor points is used to record the data. After the simulation is complete, random Gaussian noise is added using <code><a href="addNoise.html">addNoise</a></code> to give a signal to noise ratio of 40dB (based on the peak of the recorded signal and the root-mean-squared noise level).</p>

<pre class="codeinput">
<span class="comment">% run the forward simulation</span>
sensor_data = kspaceFirstOrder2D(kgrid, medium, source, sensor, input_args{:});

<span class="comment">% add noise to the recorded sensor data</span>
signal_to_noise_ratio = 40;	<span class="comment">% [dB]</span>
sensor_data = addNoise(sensor_data, signal_to_noise_ratio, 'peak');
</pre>

<a name="heading3"></a>
<h2>Selecting the regularisation parameters</h2>

<p>To compensate for the acoustic attenuation in the forward problem, the absorption parameter must be reversed in sign (the frequency content must grow rather than decay). This is achieved by setting <code>medium.alpha_sign</code> to <code>[-1, 1]</code>. The two inputs control the sign of the absorption and dispersion parameters, respectively. The dispersion parameter (which controls the dependence of the sound speed on frequency) should not be reversed because if the high frequency components of the wave field have travelled to the detector faster than the low frequency components (as is the case for photoacoustic signals in biological tissue), they again need to travel faster than the low frequency components in the time reversal reconstruction to regain their initial position within the medium.</p>

<p>If a photoacoustic image is reconstructed from real or noisy data with the absorption parameter reversed, the high frequency content (where the signal to noise ratio is typically quite low) can quickly grow to mask the low frequency content, effectively obscuring the desired features within the reconstructed image. This is because, as the waves propagate, the high frequencies are increased at a much faster rate (recall the acoustic attenuation in soft biological tissue follows a frequency power law). To avoid this, the absorption and dispersion parameters are filtered using a frequency domain Tukey window. This regularises the reconstruction (in effect, stops it from 'blowing up') by restricting the range of frequencies that are allowed to grow. The filter is applied to the absorption parameters by assigning it to <code>medium.alpha_filter</code>. The filter must have the same number of dimensions and grid points as <code>kgrid.k</code>.</p>

<p>To choose an appropriate filter cutoff frequency, the average power spectrum of the simulated sensor data is computed. This is plotted below (shown in black), along with the power spectrum of the signal recorded at the first sensor point (shown in red). In this example, the filter cutoff frequency (shown as the left dashed line in the figure) is chosen based on the noise floor observable in the power spectrum. The maximum frequency supported by the grid is also shown (right dashed line in the figure).</p>

<img vspace="5" hspace="5" src="images/example_pr_2D_tr_absorption_compensation_01.png" style="width:560px;height:420px;" alt="">

<p>Here, a symmetric Tukey window (or tapered cosine window) scaled to a particular frequency cutoff is used so that the correct power law absorption and dispersion characteristics are maintained within the filter pass band. The filter is created using <code><a href="getAlphaFilter.html">getAlphaFilter</a></code> which creates the filter using <code><a href="getWin.html">getWin</a></code> via rotation.</p> 

<pre class="codeinput">
<span class="comment">% define the cutoff frequency for the filter</span>
f_cutoff = 3e6;                     <span class="comment">% [Hz]</span>

<span class="comment">% create the filter to regularise the absorption parameters</span>
medium.alpha_filter = getAlphaFilter(kgrid_recon, medium, f_cutoff);

<span class="comment">% reverse the sign of the absorption proportionality coefficient</span>
medium.alpha_sign = [-1, 1];        <span class="comment">% [absorption, dispersion];</span>
</pre>

<p>A visualisation of the filter is given below.</p>

<img vspace="5" hspace="5" src="images/example_pr_2D_tr_absorption_compensation_02.png" style="width:560px;height:420px;" alt="">

<a name="heading4"></a>
<h2 class="title">Time reversal image reconstruction</h2>

<p>The initial pressure distribution, and the reconstructions with and without compensation for acoustic attenuation are shown below. Profiles through x = 0 are also shown for comparison. Without correction for attenuation, the edges of the reconstructed pressure become blurred, and the overall magnitude is reduced. When compensation for acoustic attenuation is included into the reconstruction, the sharpness and magnitude of the reconstructed pressure is considerably improved. Due to the band-limited frequency response of the sensors (high frequencies are not detected because their magnitudes are below the noise floor of the added noise), the reconstruction is unable to completely recover the frequency content of the initial pressure distribution. This results in Gibbs phenomenon appearing in the reconstruction where there are sharp changes in the initial pressure.</p>

<img vspace="5" hspace="5" src="images/example_pr_2D_tr_absorption_compensation_03.png" style="width:560px;height:525px;" alt="">
<img vspace="5" hspace="5" src="images/example_pr_2D_tr_absorption_compensation_04.png" style="width:560px;height:420px;" alt="">

</div></body></html>